home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_7 / a7_6.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.3 KB  |  103 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 7.6 (Gauss-Legendre Quadrature).
  9. % Section    7.5,    Gauss-Legendre Integration, Page 397
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements Gaussian quadrature.
  13.  
  14. % The function f(x) is integrated over the interval [a,b].
  15.  
  16. %    Define and store the function f(x)    in the M-file  f.m 
  17. % function y = f(x)
  18. % y = 13.*(x - x.^2).*exp(-3.*x./2);
  19.  
  20. delete f.m
  21. diary f.m; disp('function y = f(x)');...
  22.            disp('y = 13.*(x - x.^2).*exp(-3.*x./2);');...
  23. diary off;
  24.  
  25. f(0); % Remark. f.m gauss.m gauraw.m gauaw.mat are used for Algorithm 7.6
  26.  
  27. pause    % Press any key to see the graph y = f(x).
  28.  
  29. % The abscissas and weights are loaded from the file gauaw.mat
  30.  
  31. % gauraw.m contains the instructions to read this file.
  32.  
  33. clc; clg;
  34. a = 0;
  35. b = 4;
  36. h = (b-a)/200;
  37. X1 = a:h:b;
  38. Y1 = f(X1);
  39. c = -1.5;
  40. d = 2;
  41. axis([a b c d]);...
  42. plot(X1,Y1,'g');...
  43. hold on;...
  44. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  45. xlabel('x');...
  46. ylabel('y');...
  47. title('The curve y = f(x) over [a,b].');...
  48. grid;...
  49. axis;...
  50. hold off;...
  51. shg; pause    % Press any key to continue.
  52.  
  53. clc;
  54.  
  55.                 % Be patient for a moment.
  56.  
  57. [A,W] = gauraw; % Loading the abscissas and weights.
  58.  
  59.  
  60. pause    % Press any key to continue.
  61.  
  62. clc;
  63.  
  64. % Gaussian quadrature is applied over the interval [a,b].
  65.  
  66. %    Place the endpoints of [a,b] in  a  and  b.
  67.  
  68. % Place the tolerance in toler.
  69.  
  70. a  = 0;
  71. b  = 4;
  72. toler = 0.00001;
  73.  
  74. [Q,N,X,quad,err] = gauss('f',a,b,A,W,toler);
  75.  
  76. pause    % Press any key to continue.
  77.  
  78. clg;
  79. Y = f(X);
  80. Z = zeros(1,length(X));
  81. c = -1.5;
  82. d = 2;
  83. axis([a b c d]);...
  84. plot(X1,Y1,'g',X,Y,'or',X,Z,'+r');...
  85. hold on;...
  86. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  87. xlabel('x');...
  88. ylabel('y');...
  89. title('Gaussian quadrature.');...
  90. grid;...
  91. axis;...
  92. hold off;...
  93. shg; pause    % Press any key to continue.
  94.  
  95. values = [N;Q];
  96. Mx1 = 'The Gaussian quadrature approximations: ';
  97. Mx2 = '     No. of nodes       quadrature value';
  98. Mx3 = 'The approximate value of the integral is:';
  99. Mx4 = '   quadrature value   +- error bound';
  100. clc,echo off,diary output,...
  101. disp(Mx1),disp(''),disp(Mx2),disp(values'),...
  102. disp(Mx3),disp(Mx4),disp([quad err]),,diary off,echo on
  103.